# -*- coding: utf-8 -*-
"""
PBAR_ZspecCalibrationData2ndSetFirstLook.py

Created on Thu Mar 13 14:19:24 2014

@author: jkwong
"""

import numpy as np
import PBAR_Zspec, os
import matplotlib.pyplot as plt

# Set useful plot variables
plotColors = ['r', 'b', 'g', 'm', 'c', 'y', 'k'] * 10
lineStyles = ['-', '-.', ':', '_', '|'] *  10
markerTypes = ['.', 'o', 'v', '^', '<', '>', '1', '2', '3', '4', 's', 'p', '*', 'h', 'H', '+', 'x', 'D', 'd']
markerTypes = markerTypes * 2

# where all the data files are
basedir = r'C:\Users\jkwong\Documents\Work\PBAR\data4\Mar-files'
plotPath = r'C:\Users\jkwong\Documents\Work\PBAR\Plots\20140326'
# Good/bad detector masks
#badZspecList = np.array([1,2,3,4,5,6,7,8,19,20,26,62,68,76,80,121,122,123,123,125,126,127,128,129,130,131,132,133,134,135,136])
badZspecList = np.array([1,2,3,4,5,6,7,8,19,20,26,62,68,76,80,121,122,123,123,124,125,126,127,128,129,130,131,132,133,134,135,136])

badZspecMask = np.zeros(136) == 1
badZspecMask[badZspecList - 1] = True
goodZspecMask = ~badZspecMask

# read in the Zspec coordinates
zspecCoordinateFilename = 'ZSpec dimension from Source Spot 11-4-13.csv'
zspecCoordinateFullFilename = os.path.join(basedir,zspecCoordinateFilename)

temp = np.genfromtxt(zspecCoordinateFullFilename, delimiter = ',', skip_header = 1)

zspecCoordinate = {}
zspecCoordinate['x'] = temp[:,2]
zspecCoordinate['y'] = temp[:,3]
zspecCoordinate['block'] = temp[:,1]
# 
zspecMountAngle = 34. * np.pi/180.

zspecCoordinate['thetaRad'] = np.arctan(zspecCoordinate['y']/zspecCoordinate['x'] )
zspecCoordinate['thetaDeg'] = np.arctan(zspecCoordinate['y']/zspecCoordinate['x'] ) * 180 / 3.1415926535
zspecCoordinate['thetaDegOffset'] = zspecCoordinate['thetaDeg'] - 34.

# theta = 0 is horizontal
thicknessCorrection = 1 / np.cos(zspecCoordinate['thetaRad'] - zspecMountAngle)

# detector lists, descriptions, etc.
prefixList = []
#prefixList = ['ec76', 'ec77', 'ec78', 'ec79', 'ec80', 'ec81', 'ec82', 'ec83', \
#    'ec84', 'ec85', 'ec86', 'ec87', 'ec88', 'ec89', 'ec90', 'ec91', 'ec92', \
#    'ec93', 'ec94', 'ec95', 'ec96', 'ec97', 'ec98', 'ec99', 'ed00', 'ed01', 'ed02'\
#    'ed03', 'ed04', 'ed05', 'ed06']

prefix = 'ec'
for i in xrange(76,100):
    prefixList.append('%s%02d' %(prefix, i))

prefix = 'ed'
for i in xrange(0,85):
    prefixList.append('%s%02d' %(prefix, i))
# list of datasets to remove from the list
prefixRemoveList = ['ed13', 'ed17']
for p in prefixRemoveList:
    prefixList.remove(p)

# Load summary data
infoFilename = os.path.join(basedir,'DatasetSummary2ndSet.txt')

# load dataset information
(datasetDescription, datasetAcquisitionTime, \
    datasetTime, datasetTimeNum, datasetTimeStr) = \
    PBAR_Zspec.GetDatasetInformation(infoFilename, prefixList)

# pulse rate
pulserateList = np.ones(len(prefixList)) * 250.
pulserateList[0] = 60.
pulserateList[prefixList.index('ec76')] = 60.
pulserateList[prefixList.index('ec84')] = 60.
pulserateList[prefixList.index('ec85')] = 60.

# dataset groups
datasetGroups = {}
datasetGroupsWidth = {}

#datasetGroups['CC'] = ('')
# 5, 4, 6, 7, 8
datasetGroups['Pb'] = ['ec79', 'ed69', 'ec78', 'ed70', 'ec80', 'ed71',  'ec81', 'ed72', 'ec82', 'ed73']
datasetGroupsWidth['Pb'] = np.zeros( (len(datasetGroups['Pb']), 2) )
datasetGroupsWidth['Pb'][:,0] = [4,   4.5,  5.0, 5.5,  6.0,     6.5,    7.0,   7.5,    8.0,  8.5]

datasetGroups['Pb_0'] = [        'ec79', 'ec78', 'ec80', 'ec81', 'ec82']
datasetGroupsWidth['Pb_0'] = np.zeros( (len(datasetGroups['Pb_0']), 2) )
datasetGroupsWidth['Pb_0'][:,0] = [4.,     5.,    6.,     7.,     8.]

datasetGroups['Pb_1'] = [       'ed69', 'ed70', 'ed71', 'ed72',  'ed73']
datasetGroupsWidth['Pb_1'] = np.zeros( (len(datasetGroups['Pb_1']), 2) )
datasetGroupsWidth['Pb_1'][:,0] = [4.5,   5.5,    6.5,     7.5,      8.5]


# 8, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17
datasetGroups['Fe'] = ['ec87', 'ec88', 'ed61', 'ec89', 'ed62', 'ec90', 'ed63', 'ec91', 'ed64', 'ec92', 'ed65']
datasetGroupsWidth['Fe'] = np.zeros( (len(datasetGroups['Fe']), 2) )
datasetGroupsWidth['Fe'][:,0] = [8,  8,     9,       10,     11,     12,     13,     14,     15,     16,     17]

# 8, 8, 10, 12, 14, 16
datasetGroups['Fe_0'] = ['ec87', 'ec88', 'ec89', 'ec90', 'ec91', 'ec92']
datasetGroupsWidth['Fe_0'] = np.zeros( (len(datasetGroups['Fe_0']), 2) )
datasetGroupsWidth['Fe_0'][:,0] = [8,   8,     10,     12,     14,    16]

# 9, 11, 13, 15, 17
datasetGroups['Fe_1'] = ['ed61', 'ed62', 'ed63', 'ed64', 'ed65']
datasetGroupsWidth['Fe_1'] = np.zeros( (len(datasetGroups['Fe_1']), 2) )
datasetGroupsWidth['Fe_1'][:,0] = [9,  11,     13,     15,     17]


# first set: 16, 20, 24, 28  second set: 24, 26, 28, 30, 32, 34, 38
#datasetGroups['Al'] = ['ec94', 'ec95', 'ec96', 'ec97', 'ed52', 'ed53', 'ed54', 'ed55', 'ed56', 'ed57', 'ed58']
datasetGroups['Al'] = ['ec94', 'ec95', 'ec96', 'ed52', 'ed53', 'ec97', 'ed54', 'ed55', 'ed56', 'ed57', 'ed58']
datasetGroupsWidth['Al'] = np.zeros( (len(datasetGroups['Al']), 2) )
datasetGroupsWidth['Al'][:,0] = [16,  20,    24,     24,     26,    28,      28,     30,      32,    34,    38]


datasetGroups['Al_0'] = ['ec94', 'ec95', 'ec96', 'ec97']
datasetGroupsWidth['Al_0'] = np.zeros( (len(datasetGroups['Al_0']), 2) )
datasetGroupsWidth['Al_0'][:,0] = [16,  20,    24,     28]

datasetGroups['Al_1'] = ['ed52', 'ed53', 'ed54', 'ed55', 'ed56', 'ed57', 'ed58']
datasetGroupsWidth['Al_1'] = np.zeros( (len(datasetGroups['Al_1']), 2) )
datasetGroupsWidth['Al_1'][:,0] = [24, 26,     28,     30,      32,    34,    38]

# calibration and in between datasets
datasetGroups['CC'] = ['ec76', 'ec77', 'ec83', 'ec84', 'ec85', \
    'ec86', 'ec93', 'ec98', 'ec99', 'ed00', 'ed01', 'ed02', 'ed08', 'ed09', \
    'ed10', 'ed11', 'ed12', 'ed14', 'ed15', 'ed21', 'ed27', 'ed28', 'ed35',\
    'ed36', 'ed37', 'ed43', 'ed44', 'ed49', 'ed50', 'ed51', 'ed59', 'ed60', 'ed66', 'ed67', 'ed68', 'ed74', 'ed84']
datasetGroupsWidth['CC'] = np.zeros( (len(datasetGroups['CC']), 2) )
datasetGroupsWidth['CC'][:,0] = 4

# mixtures
# 3 inch Pb, Fe 4, 6, 8, 10 ,12
datasetGroups['Pb3Fe'] = ['ed03', 'ed04', 'ed05', 'ed06', 'ed07']
datasetGroupsWidth['Pb3Fe'] = np.zeros( (len(datasetGroups['Pb3Fe']), 2) )
datasetGroupsWidth['Pb3Fe'][:,0] = 3
datasetGroupsWidth['Pb3Fe'][:,1] = np.array([4, 6, 8, 10, 12])

# 4 inch Pb, Fe 2, 2, 4, 5, 6, 7, 8, 9, 10
datasetGroups['Pb4Fe'] = ['ed38', 'ed75', 'ed76', 'ed39', 'ed77', 'ed40', 'ed78', 'ed41', 'ed79', 'ed42']
datasetGroupsWidth['Pb4Fe'] = np.zeros( (len(datasetGroups['Pb4Fe']), 2) )
datasetGroupsWidth['Pb4Fe'][:,0] = 4
datasetGroupsWidth['Pb4Fe'][:,1] = np.array([2, 2, 3, 4, 5, 6, 7, 8, 9, 10])

# 4 inch Pb, Fe 2, 4, 6, 8, 10
datasetGroups['Pb4Fe_0'] = ['ed38', 'ed39', 'ed40', 'ed41', 'ed42']
datasetGroupsWidth['Pb4Fe_0'] = np.zeros( (len(datasetGroups['Pb4Fe_0']), 2) )
datasetGroupsWidth['Pb4Fe_0'][:,0] = 4
datasetGroupsWidth['Pb4Fe_0'][:,1] = np.array([2, 4, 6, 8, 10])

# 4 inch Pb, Fe 2, 3, 5, 7, 9
datasetGroups['Pb4Fe_1'] = ['ed75', 'ed76', 'ed77', 'ed78', 'ed79']
datasetGroupsWidth['Pb4Fe_1'] = np.zeros( (len(datasetGroups['Pb4Fe_1']), 2) )
datasetGroupsWidth['Pb4Fe_1'][:,0] = 4
datasetGroupsWidth['Pb4Fe_1'][:,1] = np.array([2, 3, 5, 7, 9])


# 5 inch Pb, Fe 2, 4, 6, 8
datasetGroups['Pb5Fe'] = ['ed45', 'ed46', 'ed47', 'ed48']
datasetGroupsWidth['Pb5Fe'] = np.zeros( (len(datasetGroups['Pb5Fe']), 2) )
datasetGroupsWidth['Pb5Fe'][:,0] = 5
datasetGroupsWidth['Pb5Fe'][:,1] = np.array([2, 4, 6, 8])

# 3 inch Pb, Al 12, 16, 20, 24
datasetGroups['Pb3Al'] = ['ed20', 'ed16', 'ed18', 'ed19']
datasetGroupsWidth['Pb3Al'] = np.zeros( (len(datasetGroups['Pb3Al']), 2) )
datasetGroupsWidth['Pb3Al'][:,0] = 3
datasetGroupsWidth['Pb3Al'][:,1] = np.array([12, 16, 20, 24])


# 4 inch Pb, Al 8, 12, 16, 20, 24
datasetGroups['Pb4Al'] = [    'ed22', 'ed80', 'ed81',   'ed23', 'ed82',  'ed24', 'ed83', 'ed25', 'ed26']
datasetGroupsWidth['Pb4Al'] = np.zeros( (len(datasetGroups['Pb4Al']), 2) )
datasetGroupsWidth['Pb4Al'][:,0] = 4
datasetGroupsWidth['Pb4Al'][:,1] = np.array([8., 8., 10.,  12., 14., 16., 18., 20., 24.])

# 4 inch Pb, Al 12, 16, 20, 24
datasetGroups['Pb4Al_0'] = ['ed22', 'ed23', 'ed24', 'ed25', 'ed26']
datasetGroupsWidth['Pb4Al_0'] = np.zeros( (len(datasetGroups['Pb4Al_0']), 2) )
datasetGroupsWidth['Pb4Al_0'][:,0] = 4
datasetGroupsWidth['Pb4Al_0'][:,1] = np.array([8, 12, 16, 20, 24])

# 4 inch Pb, Al 12, 16, 20, 24
datasetGroups['Pb4Al_1'] = ['ed80', 'ed81', 'ed82', 'ed83']
datasetGroupsWidth['Pb4Al_1'] = np.zeros( (len(datasetGroups['Pb4Al_1']), 2) )
datasetGroupsWidth['Pb4Al_1'][:,0] = 4
datasetGroupsWidth['Pb4Al_1'][:,1] = np.array([8, 10, 14, 18])


# 5 inch Pb, Al 4, 8, 12, 16, 20, 24
datasetGroups['Pb5Al'] = ['ed29', 'ed30', 'ed31', 'ed32', 'ed33', 'ed34']
datasetGroupsWidth['Pb5Al'] = np.zeros( (len(datasetGroups['Pb5Al']), 2) )
datasetGroupsWidth['Pb5Al'][:,0] = 5
datasetGroupsWidth['Pb5Al'][:,1] = np.array([4, 8, 12, 16, 20, 24])

# Between-Runs
datasetGroups['PbIB'] = ['ec76', 'ec83']


datasetGroups['FeIB'] = ['ec86', 'ec93', 'ed60', 'ed66'] # voltage was changed after ec86
datasetGroups['FeIB_0'] = ['ec86', 'ec93'] # voltage was changed after ec86
datasetGroups['FeIB_1'] = ['ed60', 'ed66'] # voltage was changed after ec86

datasetGroups['AlIB'] = ['ec93', 'ec98', 'ed51', 'ed59']
datasetGroups['AlIB_0'] = ['ec93', 'ec98'] # first set
datasetGroups['AlIB_1'] = ['ed51', 'ed59'] # second set, more attenuation


# mixed materials
datasetGroups['Pb3FeIB'] = ['ed02', 'ed08']

datasetGroups['Pb4FeIB'] = ['ed37', 'ed43', 'ed74', 'ed84']
datasetGroups['Pb4FeIB_0'] = ['ed37', 'ed43']
datasetGroups['Pb4FeIB_1'] = ['ed74', 'ed84']

datasetGroups['Pb5FeIB'] = ['ed44', 'ed49']

datasetGroups['Pb3AlIB'] = ['ed15', 'ed21']

datasetGroups['Pb4AlIB'] = ['ed21', 'ed27', 'ed74', 'ed84']
datasetGroups['Pb4AlIB_0'] = ['ed21', 'ed27']
datasetGroups['Pb4AlIB_1'] = ['ed74', 'ed84']

datasetGroups['Pb5AlIB'] = ['ed28', 'ed35']

datasetGroups['pure'] = datasetGroups['Pb'] + datasetGroups['Fe'] + datasetGroups['Al']
datasetGroups['PbFeMix'] =datasetGroups['Pb3Fe'] + datasetGroups['Pb4Fe'] + datasetGroups['Pb5Fe']
datasetGroups['PbAlMix'] =datasetGroups['Pb3Al'] + datasetGroups['Pb4Al'] + datasetGroups['Pb5Al']
datasetGroups['mix'] = datasetGroups['PbFeMix'] + datasetGroups['PbAlMix']


#datasetGroups['PbFe'] = ('ed03', 'ed04', 'ed05', 'ed06', 'ed07')
datasetGroups['IB'] = datasetGroups['PbIB'] + datasetGroups['FeIB'] + datasetGroups['AlIB'] + \
    datasetGroups['Pb3FeIB'] + datasetGroups['Pb4FeIB'] + datasetGroups['Pb5FeIB'] + \
    datasetGroups['Pb3AlIB'] + datasetGroups['Pb4AlIB'] + datasetGroups['Pb5AlIB']

# get unique list 
datasetGroups['IB'] = list(set(datasetGroups['IB']))
datasetGroups['IB'].sort()

# may be some redundancies to 

# all the possible dataset materials 
tempList = ['Pb', 'Fe', 'Al', 'CC', 'Pb3Fe', 'Pb4Fe', 'Pb5Fe', 'Pb3Al', 'Pb4Al', 'Pb5Al']
# create 1-to-1 arrays for material type and widths
datasetMaterial = []
datasetRawWidth = []

for j,p in enumerate(prefixList):
    print(j)
    for i,t in enumerate(tempList):
        if p in datasetGroups[t]:
            datasetMaterial.append(t)
            index = datasetGroups[t].index(p)
            datasetRawWidth.append(datasetGroupsWidth[t][index])
            break
        if i == (len(tempList)-1): # if not assigned a group then ''
            datasetMaterial.append('')
            print(p)

datasetRawWidth = np.array(datasetRawWidth)
datasetMaterial = np.array(datasetMaterial)


for i,p in enumerate(prefixList):
    print('%s, %d, %d, %s'  %(p, datasetRawWidth[i][0], datasetRawWidth[i][0], datasetDescription[i]) )

#  thickness of material seen by each detector
materialWidthArray = np.zeros((len(prefixList), 136  ))

for i, t in enumerate(datasetRawWidth):
    materialWidthArray[i,:] = t[0] * thicknessCorrection

# make mask for whether or not it is high z material
# Pb and CC(tungsten) are high z material
highZMask = np.zeros(len(datasetDescription)) == 1

for (index,d) in enumerate(datasetDescription):
    if (d.lower().find('coll') > 0  or d.lower().find('pb') > 0 ):
        highZMask[index] = True
highZMask[0] = True

mixMask = highZMask = np.zeros(len(datasetDescription)) == 1
for i in xrange(len(datasetDescription)):
    if filenameList[i] in datasetGroups['mix']:
        mixMask[i] = True


markerTypesByMaterial = {}
markerTypesByMaterial['Pb'] = 'v'
markerTypesByMaterial['Pb3Fe'] = '>'
markerTypesByMaterial['Pb4Fe'] = '<'
markerTypesByMaterial['Pb5Fe'] = '^'

markerTypesByMaterial['Pb3Al'] = 'd'
markerTypesByMaterial['Pb4Al'] = 'D'
markerTypesByMaterial['Pb5Al'] = 'h'

markerTypesByMaterial['CC'] = 'x'
markerTypesByMaterial['Fe'] = 's'
markerTypesByMaterial['Al'] = 'o'


# Get the indicies for each dataset in each group list
datasetGroupsIndices = {}
for k in datasetGroups.keys():
#    datasetGroupsIndices[k] = []
    temp = []
    for prefix in datasetGroups[k]:
#        datasetGroupsIndices[k].append(np.where(prefixList == prefix)[0][0])
        temp.append(prefixList.index(prefix))
    datasetGroupsIndices[k] = np.array(temp)

fullFilenameList = []
for p in prefixList:
    fullFilenameList.append(os.path.join(basedir,'%s.csv' %p))

# read in the data
dat = PBAR_Zspec.ReadZspec(fullFilenameList)

## Calculate the spectra stats
gainCorrectionMat = np.ones((len(dat), dat[0].shape[1]))

lowerBinThreshold = 7
baselineDatasetIndex = 0
spectrumBins = np.arange(256)
stats = PBAR_Zspec.CalculateStats(dat, datasetAcquisitionTime, spectrumBins, gainCorrectionMat, lowerBinThreshold, pulserateList)



# DONE LOADING DATA AND OTHER USEFUL STUFF



# Gain correction

# get the timestamp of the IB datasets
#datasetTimeNum
rateRange =  [0.5, 2.0]
rateLook = 1.0
binBaseline = 100
threshold = 0.05

acqTime60Hz = datasetAcquisitionTime *  pulserateList / 60.

binCalAll = np.zeros((len(dat), 136))
for i in xrange(len(dat)):
    (gainShift, score, binMean) = PBAR_Zspec.CheckGain(dat[i], acqTime60Hz[i], rateRange, rateLook, binBaseline, threshold)
    binCalAll[i,:] = binMean

# matrix for gain Correction across time


gainCorrection = np.tile(binCalAll[0,:],( len(dat), 1)  )/ binCalAll

gainExtrapolated = PBAR_Zspec.ExtrapolateGain(datasetTimeNum[datasetGroupsIndices['IB']], gainCorrection[datasetGroupsIndices['IB'],:], datasetTimeNum)






###############
###  PLOTS  ###
###############

##  Spectra  ##
# plot some spectra, same detector
detNumList = [20, 30, 40, 50, 60, 70, 80, 90, 100] # starts at 1
detNumList = [20, 40,  60, 80 , 100] # starts at 1

detNumList = [20, 40,  60, 80 , 100] # starts at 1
detNumList = [85] # starts at 1

plotIndexList = datasetGroupsIndices['CC']
#plotIndexList = [86, 16]
#plotIndexList = datasetGroupsIndices['Fe']
#
#plotIndexList = datasetGroupsIndices['Pb4Al']

#plotIndexList = np.arange(8,15)
#plotIndexList = [0, 1, 2, 3]

temp = [60., 250.]

plt.figure()
plt.grid()
for (dIndex, plotIndex) in enumerate(plotIndexList):
    for (detIndex, detNum) in enumerate(detNumList):
        plt.plot(np.arange(256), dat[plotIndex][:,detNum-1] / datasetAcquisitionTime[plotIndex], color = plotColors[dIndex %7], \
            linestyle =  lineStyles[dIndex/7], label = '%s, Det# %d, %s, %d' \
            %(prefixList[plotIndex], detNum, datasetDescription[plotIndex], pulserateList[plotIndex]))
plt.xlabel('Bin', fontsize = 16)
plt.ylabel('Count Rate (Hz)', fontsize = 16)
plt.legend()
plt.yscale('log')



# plot some spectra, same detector, counts per pulse
detNum = 101 # starts at 1

plotIndexList = datasetGroupsIndices['Pb']
plotIndexList = datasetGroupsIndices['Fe']
#plotIndexList = datasetGroupsIndices['Al']
#plotIndexList = datasetGroupsIndices['pure']
plotIndexList = datasetGroupsIndices['PbIB']
plotIndexList = datasetGroupsIndices['AlIB']
plotIndexList = datasetGroupsIndices['CC']
plotIndexList = datasetGroupsIndices['IB']

#plotIndexList = np.arange(8,15)
#plotIndexList = [0, 1, 2, 3]

temp = [60., 250.]

plt.figure()
plt.grid()
for (dIndex, plotIndex) in enumerate(plotIndexList):
    plt.plot(np.arange(256), dat[plotIndex][:,detNum-1] / datasetAcquisitionTime[plotIndex] / pulserateList[plotIndex], color = plotColors[dIndex %7], \
        linestyle =  lineStyles[dIndex/7], label = '%s, Det# %d, %s, %d, mean bin %3.3f' %(prefixList[plotIndex], detNum,datasetDescription[plotIndex], pulserateList[plotIndex], stats['binMean'][plotIndex,detNum-1]))
plt.xlabel('Bin', fontsize = 16)
plt.ylabel('Count Per Pulse', fontsize = 16)
plt.legend(loc = 1, prop = {'size':10})
plt.yscale('log')
plt.title('Detector %d' %detNum)


# plot some spectra, same detector

datasetGroupName = 'Pb3Fe'

plotIndexList = datasetGroupsIndices[datasetGroupName]

detNum = 50 # starts at 1

plt.figure()
plt.grid()
for (dIndex, d) in enumerate(dat):
    plt.plot(np.arange(256), d[:,detNum-1] / datasetAcquisitionTime[dIndex] / pulserateList[dIndex], color = plotColors[dIndex %7], \
        linestyle =  lineStyles[dIndex/7], label = '%s, Det# %d, %s' %(prefixList[dIndex], detNum,datasetDescription[dIndex]))
plt.xlabel('Bin', fontsize = 16)
plt.ylabel('Counts Per Pulse', fontsize = 16)
plt.legend()
plt.yscale('log')

# plot some spectra, same detector
detNum = 70 # starts at 1

plt.figure()
plt.grid()
for (dIndex, d) in enumerate(dat):
    plt.plot(np.arange(256), d[:,detNum-1] / datasetAcquisitionTime[dIndex], color = plotColors[dIndex %7], \
        linestyle =  lineStyles[dIndex/7], label = '%s, Det# %d, %s' %(prefixList[dIndex], detNum,datasetDescription[dIndex]))
plt.xlabel('Bin', fontsize = 16)
plt.ylabel('Count Rate (Hz)', fontsize = 16)
plt.legend()



# plot spectra, same material same thickness, different detector/angle

#thicknessRange = [5, 5.10]
#thicknessRange = [6, 6.10]
thicknessRange = [11.4,11.5]

materialName = 'Fe'
thicknessRange = [12.1,12.2]
#thicknessRange = [9.0,9.1]

#materialName = 'Al'
#thicknessRange = [10.,10.1]
#thicknessRange = [32.6, 32.8]
#
#materialName = 'Pb'
#thicknessRange = [5.0,5.1]
#thicknessRange = [5.5,5.6]
#thicknessRange = [7.0,7.10]
#thicknessRange = [8,8.2]
#thicknessRange = [8.45, 8.8]



normalized  = True

normalized  = False

thicknessRange_cut = (materialWidthArray >= thicknessRange[0]) & (materialWidthArray <= thicknessRange[1]) 
material_cut = np.zeros_like(thicknessRange_cut) == 1
material_cut[datasetGroupsIndices[materialName],:] = True
goodDetector_cut = np.zeros_like(thicknessRange_cut) == 1
goodDetector_cut[:,goodZspecMask] = True

cut = thicknessRange_cut & material_cut & goodDetector_cut

# datasets that have the
datasetIndices = np.where(cut.sum(1) > 0)[0]
maxDetectorPerDatasets = 2

angleTemp = []
statTemp = []


temp = [60., 250.]

plt.figure()
plt.grid()
ii = 0
for (dIndex, plotIndex) in enumerate(datasetIndices):
    detNumList = np.where(cut[plotIndex,:])[0] # detector number - 1
    numPos = 0
    numNeg = 0
    for (detIndex, detNum) in enumerate(detNumList):
        
        # skip if already plotted enough + angle detectors from this set
        if (numPos >= maxDetectorPerDatasets and zspecCoordinate['thetaDegOffset'][detNum] >= 0):
            continue
            print(detNum, zspecCoordinate['thetaDegOffset'][detNum])
        # skip if already plotted enough - angle detectors from this set        
        if (numNeg >= maxDetectorPerDatasets and zspecCoordinate['thetaDegOffset'][detNum] < 0):
            continue
            print(detNum, zspecCoordinate['thetaDegOffset'][detNum])

        
        if zspecCoordinate['thetaDegOffset'][detNum] < 0:
            linewidth = 3.0
        else:
            linewidth = 2.0
        
#        plt.plot(np.arange(256), dat[plotIndex][:,detNum] / datasetAcquisitionTime[plotIndex], color = plotColors[detIndex], \
#            linestyle =  lineStyles[dIndex], linewidth = linewidth, alpha = 0.5,  label = '%s, Det# %d, %s, %3.2f inch, %3.2f-deg, bin mean %3.2f' \
#            %(prefixList[plotIndex], detNum+1, datasetDescription[plotIndex], \
#            materialWidthArray[plotIndex,detNum], \
#            zspecCoordinate['thetaDegOffset'][detNum], stats['binMean'][plotIndex,detNum]))
        
        if normalized:
            plt.plot(np.arange(256), dat[plotIndex][:,detNum] / float(dat[plotIndex][:,detNum].sum()), color = plotColors[ii%7], \
                linestyle =  lineStyles[ii/7], linewidth = linewidth, alpha = 0.5,  label = '%s, Det# %d, %s, %3.2f inch, %3.2f-deg, bin mean %3.2f' \
                %(prefixList[plotIndex], detNum+1, datasetDescription[plotIndex], \
                materialWidthArray[plotIndex,detNum], \
                zspecCoordinate['thetaDegOffset'][detNum], stats['binMean'][plotIndex,detNum]))
        else:
            plt.plot(np.arange(256), dat[plotIndex][:,detNum] / datasetAcquisitionTime[plotIndex], color = plotColors[ii%7], \
                linestyle =  lineStyles[ii/7], linewidth = linewidth, alpha = 0.5,  label = '%s, Det# %d, %s, %3.2f inch, %3.2f-deg, bin mean %3.2f' \
                %(prefixList[plotIndex], detNum+1, datasetDescription[plotIndex], \
                materialWidthArray[plotIndex,detNum], \
                zspecCoordinate['thetaDegOffset'][detNum], stats['binMean'][plotIndex,detNum]))
            
            
        angleTemp.append(zspecCoordinate['thetaDegOffset'][detNum])
        statTemp.append(stats['binMean'][plotIndex,detNum])
        ii += 1
        
        if zspecCoordinate['thetaDegOffset'][detNum] >= 0:
            numPos += 1
        if zspecCoordinate['thetaDegOffset'][detNum] < 0:
            numNeg += 1
        
plt.xlabel('Bin', fontsize = 16)

if normalized:
    plt.ylabel('Counts, Normalized', fontsize = 16)
else:
    plt.ylabel('Count Rate (Hz)', fontsize = 16)
    
plt.title('%s, Thickness Range: %3.3f, %3.3f' %(materialName, thicknessRange[0], thicknessRange[1]))
plt.legend(loc = 1, prop = {'size':10})
plt.yscale('log')


##  Feature plots  ##


# plot a spectrum features vs detector number, one plot for each group
plotGroups = datasetGroupsIndices.keys()
plotGroups = ['Al', 'Al_0', 'Al_1', 'Fe']

plotGroups = ['Al', 'Fe','Pb', 'Pb4Al', 'Pb4Fe']

plotGroups = ['Al', 'Fe','Pb', 'Pb4Al', 'Pb4Fe', 'Pb3Al', 'Pb5Al', 'Pb3Fe', 'Pb5Fe']

#plotGroups = ['Pb3Al', 'Pb5Al', 'Pb3Fe', 'Pb5Fe']
#plotGroups = ['CC']

plotGroups = ['FeIB_0', 'FeIB_1', 'FeIB']

plotGroups = ['PbIB']

plotGroups = ['Fe_0', 'FeIB_0']

statName = 'binMean'
#statName = 'countPerPulse'
for k in plotGroups:
    plotIndexList = datasetGroupsIndices[k]

    plt.figure(num=None, figsize=(22, 12), dpi=80, facecolor='w', edgecolor='k')
    plt.grid()
    for (dIndex, plotIndex) in enumerate(plotIndexList):
        plt.plot(np.arange(1,137)[goodZspecMask], stats[statName][plotIndex,goodZspecMask], '.r', color = plotColors[dIndex %7], \
            linestyle =  lineStyles[dIndex/7], label = '%s, %s' %(prefixList[plotIndex],datasetDescription[plotIndex]))

    # plot the mean
    plt.plot(np.arange(1,137)[goodZspecMask], stats[statName][plotIndexList,:][:,goodZspecMask].mean(0), '.r', color = 'k', \
        linestyle =  '-', linewidth = 2, label = '%s' %('Mean'))
    
    plt.xlabel('Detector Number', fontsize = 20)
    plt.ylabel(statName, fontsize = 20)
    plt.title(k, fontsize = 20)
    plt.legend(loc = 1, prop = {'size':8})
    
    if savePlots:
        filename = os.path.join(plotPath, '%s_%sVsDetectorNumber.png' %(statName, k))
        plt.savefig(filename)



# plot a spectrum features vs detector number, multiple groups for each plot
plotGroups = datasetGroupsIndices.keys()
plotGroups = ['Al', 'Al_0', 'Al_1', 'Fe']
plotGroups = ['Al', 'Fe','Pb', 'Pb4Al', 'Pb4Fe']
plotGroups = ['Al', 'Fe','Pb', 'Pb4Al', 'Pb4Fe', 'Pb3Al', 'Pb5Al', 'Pb3Fe', 'Pb5Fe']
#plotGroups = ['Pb3Al', 'Pb5Al', 'Pb3Fe', 'Pb5Fe']
#plotGroups = ['CC']
plotGroups = ['FeIB_0', 'FeIB_1', 'FeIB']
plotGroups = ['PbIB']
plotGroups = ['Fe_0', 'FeIB_0']

statName = 'binMean'
#statName = 'countPerPulse'
plt.figure(num=None, figsize=(22, 12), dpi=80, facecolor='w', edgecolor='k')
plt.grid()

plotIndexList = []
for k in plotGroups:
    plotIndexList = plotIndexList + list(datasetGroupsIndices[k])
for (dIndex, plotIndex) in enumerate(plotIndexList):
    plt.plot(np.arange(1,137)[goodZspecMask], stats[statName][plotIndex,goodZspecMask], '.r', color = plotColors[dIndex %7], \
        linestyle =  lineStyles[dIndex/7], label = '%s, %s' %(prefixList[plotIndex],datasetDescription[plotIndex]))

plt.xlabel('Detector Number', fontsize = 20)
plt.ylabel(statName, fontsize = 20)
plt.title(k, fontsize = 20)
plt.legend(loc = 1, prop = {'size':8})

if savePlots:
    filename = os.path.join(plotPath, '%s_%sVsDetectorNumber.png' %(statName, k))
    plt.savefig(filename)


# plot of the std/mean of the CC, IB datasets
#plotGroups = ['PbIB']

plotGroups = ['IB']

statName = 'binMean'
#statName = 'countPerPulse'
for k in plotGroups:
    plotIndexList = datasetGroupsIndices[k]

    plt.figure(num=None, figsize=(22, 12), dpi=80, facecolor='w', edgecolor='k')
    plt.grid()

    # plot the mean
    tempMean = stats[statName][plotIndexList,:][:,goodZspecMask].mean(0)
    tempStd = stats[statName][plotIndexList,:][:,goodZspecMask].std(0)
    
    plt.plot(np.arange(1,137)[goodZspecMask], tempStd/tempMean, '.r', color = 'k', \
        linestyle =  '-', linewidth = 2, label = '%s' %('Mean'))
    
    plt.xlabel('Detector Number', fontsize = 20)
    plt.ylabel('std/mean, mean bin', fontsize = 20)
    plt.title(k, fontsize = 20)
    plt.legend(loc = 1, prop = {'size':8})
    
    



# Plot of Feature vs Angle
plotGroups = datasetGroupsIndices.keys()
plotGroups = ['Al', 'Fe','Pb', 'Pb4Al', 'Pb4Fe']
plotGroups = ['Pb3Al', 'Pb5Al', 'Pb3Fe', 'Pb5Fe']
plotGroups = ['Al', 'Fe','Pb', 'Pb4Al', 'Pb4Fe', 'Pb3Al', 'Pb5Al', 'Pb3Fe', 'Pb5Fe']
#plotGroups = ['CC']

statName = 'binMean'
statName = 'countPerPulse'
savePlots = 1
for k in plotGroups:
    plotIndexList = datasetGroupsIndices[k]
    plt.figure(num=None, figsize=(22, 12), dpi=80, facecolor='w', edgecolor='k')
    plt.grid()
    for (dIndex, plotIndex) in enumerate(plotIndexList):
        plt.plot(zspecCoordinate['thetaDegOffset'][goodZspecMask], stats[statName][plotIndex,goodZspecMask], '.r', color = plotColors[dIndex %7], \
            linestyle =  lineStyles[dIndex/7], label = '%s, %s' %(prefixList[plotIndex],datasetDescription[plotIndex]))
    plt.xlabel('Angle (degree)', fontsize = 20)
    plt.ylabel(statName, fontsize = 20)
    plt.title(k, fontsize = 20)
    plt.legend(loc = 1, prop = {'size':12})
    if savePlots:
        filename = os.path.join(plotPath, '%s_%sVsAngle.png' %(statName, k))
        plt.savefig(filename)


# Plot of feature vs abs (angle - 45)

for k in datasetGroupsIndices.keys():
    plotIndexList = datasetGroupsIndices[k]
    statName = 'binMean'
    plt.figure(num=None, figsize=(22, 12), dpi=80, facecolor='w', edgecolor='k')
    plt.grid()
    for (dIndex, plotIndex) in enumerate(plotIndexList):
        plt.plot(np.abs(zspecCoordinate['thetaDeg'][goodZspecMask] - 45), stats[statName][plotIndex,goodZspecMask], '.r', color = plotColors[dIndex %7], \
            linestyle =  lineStyles[dIndex/7], label = '%s, %s' %(prefixList[plotIndex],datasetDescription[plotIndex]))
    plt.xlabel('Angle (degree)', fontsize = 20)
    plt.ylabel(statName, fontsize = 20)
    plt.title(k, fontsize = 20)
    plt.legend(loc = 2, prop = {'size':12})


# plot mean bin vs detector number, CC, also plot mean of all the CC datasets

datasetGroupName = 'CC'
plotIndexList = datasetGroupsIndices[datasetGroupName]

statName = 'binMean'
plt.figure()
plt.grid()

tempArray = np.zeros((len(plotIndexList), 136))
for (dIndex, plotIndex) in enumerate(plotIndexList):
    tempArray[dIndex,:] = stats[statName][plotIndex,:]
    plt.plot(np.arange(1,137)[goodZspecMask], stats[statName][plotIndex,goodZspecMask], '.r', color = plotColors[dIndex %7], \
        linestyle =  lineStyles[dIndex/7], label = '%s, %s' %(prefixList[plotIndex],datasetDescription[plotIndex]))

# plot the average
plt.plot(np.arange(1,137)[goodZspecMask], tempArray[:,goodZspecMask].mean(0), '-r', linewidth = 2, label = 'Average')
        
plt.xlabel('Detector Number', fontsize = 16)
plt.ylabel(statName, fontsize = 16)
plt.legend(loc = 2, prop = {'size':8})



# plot std(binMean) / mean(binMean) vs detector number
plotIndexList = datasetGroupsIndices['CC']
statName = 'binMean'
plt.figure()
plt.grid()
temp = stats[statName][plotIndexList,:][:,goodZspecMask].std(0) / stats[statName][plotIndexList,:][:,goodZspecMask].mean(0)
plt.plot(np.arange(1,137)[goodZspecMask], temp, '.r')
plt.xlabel('Detector Number', fontsize = 16)
plt.ylabel('std(binMean)/mean(binMean)', fontsize = 16)
plt.legend(loc = 2, prop = {'size':8})


# discrimination plots
# plot of the bin mean vs count/pulse, several detectors

statName = 'binMean'
statName = 'binSkew'
detectorMask = np.zeros(136) == 1
detectorMask[50:55] = True
plt.figure()
plt.grid()

statName1 = 'countPerPulse'
statName2 = 'binSkew'
for i in xrange(len(dat)):
    cut = goodZspecMask & detectorMask
    x = stats[statName1][i,cut]
    y = stats[statName2][i,cut]
    if highZMask[i]: # high z
        markerSize = 15
        markerAlpha = 1
    else: # not high z
        markerSize = 15
        markerAlpha = 0.5
    markerType = markerTypesByMaterial[datasetMaterial[i]]
#    plt.plot(x, y, '.k', color = plotColors[i % 7] , marker = markerTypes[i/7], markersize = markerSize, \
#        alpha = markerAlpha, label = '%s, %s' %(prefixList[i], datasetDescription[i]))
    plt.plot(x, y, '.k', color = plotColors[i % 7] , marker = markerType  , markersize = markerSize, \
        alpha = markerAlpha, label = '%s, %s' %(filenameList[i], datasetDescription[i]))
        
plt.xlabel(statName1, fontsize = 16)
plt.ylabel(statName2, fontsize = 16)
plt.legend(loc = 2, prop = {'size':8})
plt.xscale('log')



# discrimination plots
# plot of the bin mean vs count/pulse, several detectors
# high z = blue
# low z = red
# PbFe = light
# PbAl = lighter

groupList = ['Pb', 'Fe', 'Al', 'Pb3Fe', 'Pb4Fe', 'Pb5Fe', 'Pb3Al', 'Pb4Al', 'Pb5Al']
#groupList = ['Pb', 'Fe', 'Al', 'Pb3Fe', 'Pb4Fe', 'Pb5Fe']

#groupList = ['Pb', 'Fe', 'Al']

plotIndexList = []
for g in groupList:
    plotIndexList = plotIndexList + list(datasetGroupsIndices[g])

detectorMask = np.zeros(136) == 1
#detectorMask[50:55] = True
#detectorMask[20:120] = True
#detectorMask[55:60] = True
#detectorMask[80:85] = True
#detectorMask[85:86] = True

#detectorMask[40:45] = True

#detectorMask[20:25] = True
detectorMask[70:71] = True

#detectorMask[21:22] = True

#plt.figure()
plt.figure(num=None, figsize=(22, 12), dpi=80, facecolor='w', edgecolor='k')
plt.grid()

statName1 = 'countPerPulse'
statName2 = 'binSkew'
#statName2 = 'binSTD'

for i in plotIndexList:
    cut = goodZspecMask & detectorMask
    x = stats[statName1][i,cut]
#    y = stats[statName2][i,cut] * gainExtrapolated[i,cut]
    y = stats[statName2][i,cut]

#    if highZMask[i]: # high z
#        markerSize = 8
#        color = 'b'
#    else: # not high z
#        markerSize = 8
#        color = 'r'
        
    # determine the color
    if filenameList[i] in datasetGroups['Fe']:
        color = 'r'
    elif filenameList[i] in datasetGroups['Al']:
        color = 'm'
    elif filenameList[i] in datasetGroups['Pb']:
        color = 'b'
    elif filenameList[i] in datasetGroups['PbFeMix']:
        color = 'b'
    elif filenameList[i] in datasetGroups['PbAlMix']:
        color = 'b'
        
    # marker alpha        
    if filenameList[i] in datasetGroups['PbFeMix']:
        markerAlpha = 0.4
    elif filenameList[i] in datasetGroups['PbAlMix']:
        markerAlpha = 0.1
    else:
        markerAlpha = 1.0
        
    markerType = markerTypesByMaterial[datasetMaterial[i]]
#    plt.plot(x, y, '.k', color = plotColors[i % 7] , marker = markerTypes[i/7], markersize = markerSize, \
#        alpha = markerAlpha, label = '%s, %s' %(prefixList[i], datasetDescription[i]))
    plt.plot(x, y, '.k', color = color , marker = markerType  , markersize = markerSize, \
        alpha = markerAlpha, label = '%s, %s' %(filenameList[i], datasetDescription[i]))
    plt.text(np.mean(x), np.mean(y), filenameList[i])
        
plt.xlabel(statName1, fontsize = 16)
plt.ylabel(statName2, fontsize = 16)
plt.legend(loc = 2, prop = {'size':8})
plt.xscale('log')

plt.title('uncorrected')
#plt.title('Corrected')
plt.ylim((0, 50))




## POLYNOMIAL FITS PLOTS

# Plot of the sum error of MEAN best fit as a function of polynomial order

plt.figure()
plt.grid()
#for varIndex, varName in enumerate(['mean', 'sigma']):
varName = 'mean'

for groupNameIndex, groupName in enumerate(groupNamesExportList):  # 'PbNOT', 'PbALL'
    for subsetNameIndex, subsetName in enumerate(['all', 'tail', 'extend']):
        plt.plot(polyOrderList, pfitError[varName][groupName][subsetName], '-k', markersize = 15, \
            marker = markerTypes[subsetNameIndex], color = plotColors[groupNameIndex], \
            label = '%s, %s, %s' %(varName, groupName, subsetName))

plt.yscale('log')
plt.xlabel('Polynomial Order', fontsize = 14)
plt.ylabel('Sum Squared Error', fontsize = 14)
plt.legend()
plt.title('')


# Plot of the sum error of SIGMA best fit as a function of polynomial order
plt.figure()
plt.grid()

#for varIndex, varName in enumerate(['mean', 'sigma']):
varName = 'sigma'

for groupNameIndex, groupName in enumerate(groupNamesExportList):  # 'PbNOT', 'PbALL'
    for subsetNameIndex, subsetName in enumerate(['all', 'tail', 'extend']):
        plt.plot(polyOrderList, pfitError[varName][groupName][subsetName], '-k', markersize = 15, \
            marker = markerTypes[subsetNameIndex], color = plotColors[groupNameIndex], \
            label = '%s, %s, %s' %(varName, groupName, subsetName))

plt.yscale('log')
plt.xlabel('Polynomial Order', fontsize = 14)
plt.ylabel('Sum Squared Error', fontsize = 14)
plt.legend()
plt.title('')




plt.figure()
plt.grid()
#for varIndex, varName in enumerate(['mean', 'sigma']):
varName = 'sigma'

for groupNameIndex, groupName in enumerate(groupNamesExportList):  # 'PbNOT', 'PbALL'
    for subsetNameIndex, subsetName in enumerate(['all', 'tail', 'extend']):
        plt.plot(polyOrderList, pfitError[varName][groupName][subsetName], '-k', markersize = 15, \
            marker = markerTypes[subsetNameIndex], color = plotColors[groupNameIndex], \
            label = '%s, %s, %s' %(varName, groupName, subsetName))

plt.yscale('log')
plt.xlabel('Polynomial Order', fontsize = 14)
plt.ylabel('Sum Squared Error', fontsize = 14)
plt.legend()
plt.title('')


plt.figure()
plt.grid()
for varIndex, varName in enumerate(['mean', 'sigma']):
    for groupNameIndex, groupName in enumerate(groupNamesExportList):  # 'PbNOT', 'PbALL'
        for subsetNameIndex, subsetName in enumerate(['all', 'tail', 'extend']):
            plt.plot(polyOrderList, pfitError[varName][groupName][subsetName], '.k', fontsize = 15, \
                marker = markerTypes[varIndex], color = plotColors[groupNameIndex], \
                linestyle = lineStyles[subsetNameIndex], label = '%s, %s, %s' %(varName, groupName, subsetName))
plt.xlabel('Polynomial Order', fontsize = 14)
plt.ylabel('Sum Squared Error', fontsize = 14)
plt.legend()
plt.title('')




# Mean, Sigma Fit Plots


# sigma

polyOrderIndex = 9
polyOrderTailIndex = 4
polyOrderExtendIndex = 13
varName = 'sigma'

plt.figure()
plt.grid()

groupName = 'PbALL'
subsetName = 'extend'
plt.plot(fitCountBinCentersList, fitParams[groupName][:,2], '-xg', markersize = 10, linewidth = 2, markeredgewidth = 2)
x = np.arange(0, 400)
plt.plot(x, np.polyval(pfit[varName][groupName]['all'][polyOrderIndex], x), '--g', \
    linewidth = 3, label = '%s, %s, All points Order %d, %3.3f' %(groupName, varName, polyOrderIndex, pfitError[varName][groupName]['all'][polyOrderIndex]))
plt.plot(x, np.polyval(pfit[varName][groupName]['tail'][polyOrderTailIndex], x), '-g', \
    linewidth = 3, label = '%s, %s, Tail points, Order %d, %3.3f' %(groupName, varName, polyOrderTailIndex, pfitError[varName][groupName]['tail'][polyOrderTailIndex]))
plt.plot(x, np.polyval(pfit[varName][groupName]['extend'][polyOrderExtendIndex], x), '-.g', \
    linewidth = 3, label = '%s, %s, Include extra tail points, Order %d, %3.3f' %(groupName, varName, polyOrderExtendIndex, pfitError[varName][groupName]['extend'][polyOrderExtendIndex]))

groupName = 'PbNOT'
subsetName = 'extend'
plt.plot(fitCountBinCentersList, fitParams[groupName][:,2], '-xr', markersize = 10, linewidth = 2, markeredgewidth = 2)
x = np.arange(0, 400)
plt.plot(x, np.polyval(pfit[varName][groupName]['all'][polyOrderIndex], x), '--r', \
    linewidth = 3, label = '%s, %s, All points Order %d, %3.3f' %(groupName, varName, polyOrderIndex, pfitError[varName][groupName]['all'][polyOrderIndex]))
plt.plot(x, np.polyval(pfit[varName][groupName]['tail'][polyOrderTailIndex], x), '-r', \
    linewidth = 3, label = '%s, %s, Tail points, Order %d, %3.3f' %(groupName, varName, polyOrderTailIndex, pfitError[varName][groupName]['tail'][polyOrderTailIndex]))
plt.plot(x, np.polyval(pfit[varName][groupName]['extend'][polyOrderExtendIndex], x), '-.r', \
    linewidth = 3, label = '%s, %s, Include extra tail points, Order %d, %3.3f' %(groupName, varName, polyOrderExtendIndex, pfitError[varName][groupName]['extend'][polyOrderExtendIndex]))

plt.ylim((0, 5))
plt.legend()
plt.xlabel('Count')
plt.ylabel('Distance %s' % varName)
plt.title('')



# Mean

polyOrderIndex = 9
polyOrderTailIndex = 4
polyOrderExtendIndex = 22
varName = 'mean'


plt.figure()
plt.grid()

groupName = 'PbALL'
subsetName = 'extend'
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1], '-xg', markersize = 10, linewidth = 2, markeredgewidth = 2)
x = np.arange(0, 400)
plt.plot(x, np.polyval(pfit[varName][groupName]['all'][polyOrderIndex], x), '--g', \
    linewidth = 3, label = '%s, %s, All points Order %d, %3.3f' %(groupName, varName, polyOrderIndex, pfitError[varName][groupName]['all'][polyOrderIndex]))
plt.plot(x, np.polyval(pfit[varName][groupName]['tail'][polyOrderTailIndex], x), '-g', \
    linewidth = 3, label = '%s, %s, Tail points, Order %d, %3.3f' %(groupName, varName, polyOrderTailIndex, pfitError[varName][groupName]['tail'][polyOrderTailIndex]))
plt.plot(x, np.polyval(pfit[varName][groupName]['extend'][polyOrderExtendIndex], x), '-.g', \
    linewidth = 3, label = '%s, %s, Include extra tail points, Order %d, %3.3f' %(groupName, varName, polyOrderExtendIndex, pfitError[varName][groupName]['extend'][polyOrderExtendIndex]))

groupName = 'PbNOT'
subsetName = 'extend'
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1], '-xr', markersize = 10, linewidth = 2, markeredgewidth = 2)
x = np.arange(0, 400)
plt.plot(x, np.polyval(pfit[varName][groupName]['all'][polyOrderIndex], x), '--r', \
    linewidth = 3, label = '%s, %s, All points Order %d, %3.3f' %(groupName, varName, polyOrderIndex, pfitError[varName][groupName]['all'][polyOrderIndex]))
plt.plot(x, np.polyval(pfit[varName][groupName]['tail'][polyOrderTailIndex], x), '-r', \
    linewidth = 3, label = '%s, %s, Tail points, Order %d, %3.3f' %(groupName, varName, polyOrderTailIndex, pfitError[varName][groupName]['tail'][polyOrderTailIndex]))
plt.plot(x, np.polyval(pfit[varName][groupName]['extend'][polyOrderExtendIndex], x), '-.r', \
    linewidth = 3, label = '%s, %s, Include extra tail points, Order %d, %3.3f' %(groupName, varName, polyOrderExtendIndex, pfitError[varName][groupName]['extend'][polyOrderExtendIndex]))

plt.ylim((-10, 20))
plt.legend()
plt.xlabel('Count')
plt.ylabel('Distance %s' % varName)
plt.title('')


# SCATTER PLOT OF DIST VS COUNT AND THE BEST FIT LINE using best w for data points of all counts
# used for monthly report 4/12/2014

polyOrderIndex = 9
polyOrderTailIndex = 4
polyMeanOrderExtendIndex = 22
polySigmaOrderExtendIndex = 13

sigmaFraction = 2.5

plt.figure()
plt.grid()
for i in xrange(2):
    cutt = y_All == i
    plt.plot(count_All[cutt], dist0Best[cutt], '.k', color = plotColors[i], alpha = 0.2, markersize = 10, label = groupNamesExportList[i])

groupName = 'PbNOT'
groupName = 'PbALL'

varName = 'mean'
subsetName = 'extend'
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1], '-xb', markersize = 10, linewidth = 1, markeredgewidth = 2)
# sigma
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] + sigmaFraction*fitParams[groupName][:,2], '-xb', markersize = 10, linewidth = 2, markeredgewidth = 2)
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] - sigmaFraction*fitParams[groupName][:,2], '-xb', markersize = 10, linewidth = 2, markeredgewidth = 2)


# plot best fit lines
x = np.arange(0, 400)
# mean
plt.plot(x, np.polyval(pfit[varName][groupName]['extend'][polyMeanOrderExtendIndex], x), '-.b', \
    linewidth = 3, label = '%s, %s, Include extra tail points, Order %d, Error %3.3f' %(groupName, varName, polyMeanOrderExtendIndex, pfitError[varName][groupName]['extend'][polyMeanOrderExtendIndex]))
# sigma
plt.plot(x, np.polyval(pfit['mean'][groupName]['extend'][polyMeanOrderExtendIndex], x) + sigmaFraction*np.polyval(pfit['sigma'][groupName]['extend'][polySigmaOrderExtendIndex], x), '-.b', \
    linewidth = 3)
plt.plot(x, np.polyval(pfit['mean'][groupName]['extend'][polyMeanOrderExtendIndex], x) - sigmaFraction*np.polyval(pfit['sigma'][groupName]['extend'][polySigmaOrderExtendIndex], x), '-.b', \
    linewidth = 3)


groupName = 'PbNOT'

plt.plot(fitCountBinCentersList, fitParams[groupName][:,1], '-xr', markersize = 10, linewidth = 1, markeredgewidth = 2)
# sigma
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] + sigmaFraction*fitParams[groupName][:,2], '-xr', markersize = 10, linewidth = 2, markeredgewidth = 2)
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] - sigmaFraction*fitParams[groupName][:,2], '-xr', markersize = 10, linewidth = 2, markeredgewidth = 2)

#xPolygon = np.hstack((fitCountBinCentersList, np.flipud(fitCountBinCentersList)))
#
#yPolygon = np.hstack((fitParams[groupName][:,1] + 2.5*fitParams[groupName][:,2], np.flipud(fitParams[groupName][:,1] - 2.5*fitParams[groupName][:,2])))
#
#polygon = Polygon(np.random.rand(N,2), True)
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] + sigmaFraction*fitParams[groupName][:,2], '-xr', markersize = 10, linewidth = 2, markeredgewidth = 2)
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] - sigmaFraction*fitParams[groupName][:,2], '-xr', markersize = 10, linewidth = 2, markeredgewidth = 2)

# plot best fit lines
x = np.arange(0, 400)
# mean
plt.plot(x, np.polyval(pfit[varName][groupName]['extend'][polyMeanOrderExtendIndex], x), '-.r', \
    linewidth = 3, label = '%s, %s, Include extra tail points, Order %d, Error %3.3f' %(groupName, varName, polyMeanOrderExtendIndex, pfitError[varName][groupName]['extend'][polyMeanOrderExtendIndex]))
# sigma
plt.plot(x, np.polyval(pfit['mean'][groupName]['extend'][polyMeanOrderExtendIndex], x) + sigmaFraction*np.polyval(pfit['sigma'][groupName]['extend'][polySigmaOrderExtendIndex], x), '-.r', \
    linewidth = 3)
plt.plot(x, np.polyval(pfit['mean'][groupName]['extend'][polyMeanOrderExtendIndex], x) - sigmaFraction*np.polyval(pfit['sigma'][groupName]['extend'][polySigmaOrderExtendIndex], x), '-.r', \
    linewidth = 3)

plt.ylim((-20, 20))
plt.legend()
plt.xlabel('Count')
plt.ylabel('Distance')
plt.title('')




# SCATTER PLOT OF DIST VS COUNT AND THE BEST FIT LINE using best w for data points of all counts
# PB offset removed
# used for monthly report 4/12/2014

polyOrderIndex = 9
polyOrderTailIndex = 4
polyMeanOrderExtendIndex = 22
polySigmaOrderExtendIndex = 13

sigmaFraction = 2.5

plt.figure()
plt.grid()
for i in xrange(2):
    cutt = y_All == i
    plt.plot(count_All[cutt], dist0Best[cutt] - np.polyval(pfit['mean']['PbALL']['extend'][polyMeanOrderExtendIndex], count_All[cutt]), '.k', \
        color = plotColors[i], alpha = 0.2, markersize = 10, label = groupNamesExportList[i])

#groupName = 'PbNOT_0'
groupName = 'PbALL'

varName = 'mean'
subsetName = 'extend'

distOffset = np.polyval(pfit['mean']['PbALL']['extend'][polyMeanOrderExtendIndex], fitCountBinCentersList)

plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] - distOffset, '-xb', \
    markersize = 10, linewidth = 1, markeredgewidth = 2)
# sigma
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] + sigmaFraction*fitParams[groupName][:,2] - distOffset, '-xb', \
    markersize = 10, linewidth = 2, markeredgewidth = 2)
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] - sigmaFraction*fitParams[groupName][:,2] - distOffset, '-xb', \
    markersize = 10, linewidth = 2, markeredgewidth = 2)

# plot best fit lines
x = np.arange(0, 400)
# mean
plt.plot(x, np.polyval(pfit[varName][groupName]['extend'][polyMeanOrderExtendIndex], x), '-.b', \
    linewidth = 3, label = '%s, %s, Include extra tail points, Order %d, Error %3.3f' %(groupName, varName, polyMeanOrderExtendIndex, pfitError[varName][groupName]['extend'][polyMeanOrderExtendIndex]), alpha = 0.5)
# sigma
plt.plot(x, +sigmaFraction*np.polyval(pfit['sigma'][groupName]['extend'][polySigmaOrderExtendIndex], x), '-.b', \
    linewidth = 3, alpha = 0.5)
plt.plot(x, -sigmaFraction*np.polyval(pfit['sigma'][groupName]['extend'][polySigmaOrderExtendIndex], x), '-.b', \
    linewidth = 3, alpha = 0.5)


groupName = 'PbNOT'

distOffset = np.polyval(pfit['mean']['PbALL']['extend'][polyMeanOrderExtendIndex], fitCountBinCentersList)

plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] - distOffset, '-xr', \
    markersize = 10, linewidth = 1, markeredgewidth = 2)
# sigma
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] + sigmaFraction*fitParams[groupName][:,2] - distOffset, '-xr', \
    markersize = 10, linewidth = 2, markeredgewidth = 2)
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] - sigmaFraction*fitParams[groupName][:,2] - distOffset, '-xr', \
    markersize = 10, linewidth = 2, markeredgewidth = 2)

#polygon = Polygon(np.random.rand(N,2), True)
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] + sigmaFraction*fitParams[groupName][:,2] - distOffset, '-xr', markersize = 10, linewidth = 2, markeredgewidth = 2, alpha = 0.5)
plt.plot(fitCountBinCentersList, fitParams[groupName][:,1] - sigmaFraction*fitParams[groupName][:,2] - distOffset, '-xr', markersize = 10, linewidth = 2, markeredgewidth = 2, alpha = 0.5)


# plot best fit lines
x = np.arange(0, 400)
distOffset = np.polyval(pfit['mean']['PbALL']['extend'][polyMeanOrderExtendIndex], x)

# mean
plt.plot(x, np.polyval(pfit[varName][groupName]['extend'][polyMeanOrderExtendIndex], x) - distOffset, '-.r', \
    linewidth = 3, label = '%s, %s, Include extra tail points, Order %d, Error %3.3f' %(groupName, varName, polyMeanOrderExtendIndex, pfitError[varName][groupName]['extend'][polyMeanOrderExtendIndex]), alpha = 0.5)
# sigma
plt.plot(x, np.polyval(pfit['mean'][groupName]['extend'][polyMeanOrderExtendIndex], x) + sigmaFraction*np.polyval(pfit['sigma'][groupName]['extend'][polySigmaOrderExtendIndex], x) - distOffset, '-.r', \
    linewidth = 3, alpha = 0.5)
plt.plot(x, np.polyval(pfit['mean'][groupName]['extend'][polyMeanOrderExtendIndex], x) - sigmaFraction*np.polyval(pfit['sigma'][groupName]['extend'][polySigmaOrderExtendIndex], x) - distOffset, '-.r', \
    linewidth = 3, alpha = 0.5)

plt.ylim((-15, 15))
plt.legend()
plt.xlabel('Count')
plt.ylabel('Distance')
plt.title('')


# SCATTER PLOT OF DIST VS COUNT AND THE BEST FIT LINE using best w for data points of all counts
# PB offset remove, no lines

polyOrderIndex = 9
polyOrderTailIndex = 4
polyMeanOrderExtendIndex = 22
polySigmaOrderExtendIndex = 13

sigmaFraction = 2.5

plt.figure()
plt.grid()
for i in xrange(2):
    cutt = y_All == i
    plt.plot(count_All[cutt], dist0Best[cutt] - np.polyval(pfit['mean']['PbALL']['extend'][polyMeanOrderExtendIndex], count_All[cutt]), '.k', \
        color = plotColors[i], alpha = 0.2, markersize = 10, label = groupNamesExportList[i])

plt.ylim((-15, 15))
plt.legend()
plt.xlabel('Count')
plt.ylabel('Distance')
plt.title('')








#############################
###  MISCELLANEOUS PLOTS  ###
#############################


# Plot of detector positions

detectorAt45 = np.argmin( abs(zspecCoordinate['thetaRad'] - zspecMountAngle))

plt.figure()
plt.grid()
plt.plot(zspecCoordinate['x'], zspecCoordinate['y'], '.k')
plt.plot(zspecCoordinate['x'][detectorAt45], zspecCoordinate['y'][detectorAt45], '.r')

plt.plot([0,275*np.cos(zspecMountAngle)], [0,275*np.sin(zspecMountAngle)], '-r')
plt.plot(0,0, 'xk')

plt.xlabel('x', fontsize = 16)
plt.ylabel('y', fontsize = 16)

plt.title('Zspec Detector Locations')
plt.axis((-15, 225, -15, 225))


# Plot of thickness correction factor vs thickness
plt.figure()
plt.grid()
plt.plot(np.arange(136)+1, thicknessCorrection, '.-k')
plt.xlabel('Detector Number', fontsize = 16)
plt.ylabel('Thickness Correction Factor', fontsize = 16)
plt.title('Thickness Correction Factor vs Detector Number', fontsize = 16)


##  thickness vs detector

materialName = 'Pb'

datasetIndices = datasetGroupsIndices[materialName]

plt.figure()
plt.grid()
for (dIndex, plotIndex) in enumerate(datasetIndices):    
    plt.plot(np.arange(136)+1,  materialWidthArray[plotIndex,:], '.k', color = plotColors[dIndex], \
        label = '%s, %s, %3.2f' \
        %(prefixList[plotIndex], datasetDescription[plotIndex], \
        datasetRawWidth[plotIndex,0]))

for (dIndex, plotIndex) in enumerate(datasetIndices):
    t = datasetRawWidth[plotIndex,0]
    plt.plot([0, 140], [t, t], '-k')

plt.legend(loc = 1, prop = {'size':10})
plt.xlabel('Detector Number', fontsize = 16)
plt.ylabel('Thickness', fontsize = 16)
plt.title('%s, Thickness vs Detector Number' %materialName, fontsize = 16)

# histogram of thickness of materials

materialName = 'Fe'
#materialName = 'Pb'
#materialName = 'Al'

binedges = np.arange(4, 34, 0.10)
bincenters = (binedges[1:] + binedges[:-1]) / 2.
thicknessMatrix = materialWidthArray[datasetGroupsIndices[materialName],:][:,goodZspecMask]
thicknessDatasetCount = np.zeros(len(binedges) - 1)
for i in xrange(thicknessMatrix.shape[0]):
    temp = plt.hist(thicknessMatrix[i,:], bins = binedges)
    thicknessDatasetCount = thicknessDatasetCount + (temp[0] > 0)



plt.figure()
plt.grid()
plt.plot(bincenters, thicknessDatasetCount)
plt.title(materialName)
plt.xlabel('Thickness', fontsize = 16)
plt.ylabel('Dataset Count', fontsize = 16)

# plot correlation between materials datasets binmean and associated CC datasets
statName = 'binMean'
plt.figure()
plt.grid()


for index,plotIndex2 in enumerate(datasetGroupsIndices['Fe']):
    if (plotIndex >= datasetGroupsIndices['FeIB'][0] and plotIndex <= datasetGroupsIndices['FeIB'][1]):
        plotIndex1 = datasetGroupsIndices['FeIB'][0]
    else:
        plotIndex1 = datasetGroupsIndices['FeIB'][2]

    plt.plot(stats[statName][plotIndex1,:][:,goodZspecMask], stats[statName][plotIndex2,:][:,goodZspecMask], '.k',\
        color = plotColors[index%7], marker = markerTypes[index/7], label = '%s, %s' %(prefixList[plotIndex2], datasetDescription[plotIndex2]))

plt.legend()
plt.xlabel('Bin Mean CC', fontsize = 16)
plt.ylabel('Bin Mean Fe', fontsize = 16)
